//#include <vector>
//#include <iostream>
//#include <fstream>
//#include "parameters.h"
//#include "constants.h"
//#include "ljforce.h"
//#include "md.h"
//#include "pimd.h"
//#include "lscivr.h"
//#include "analysis.h"

//* main routine to run dynamics
// */
//int main2(int argc, char *argv[])
//{
//    // check input arguments @TODO
//    for (int i = 0; i < argc; i++)
//    {
//        std::cout << "argument[" << i << "] is: " << argv[i] << std::endl;
//    }
//    // @TODO write parameters from GUI

//    // set parameters for simulation, molecular system, force fields
//    std::string parameter_file = "parameters.txt";
//    std::ifstream ifs(parameter_file);
//    if (ifs.good())
//    {
//        ReadParameters(parameter_file);
//    }

//    // construct object of molecular system
//    MolSys mol;
//    // set force field
//    Force *pForce = nullptr;
//    if (gForceType == ForceTypes::kLennardJones)
//    {
//        LjForce lj(gLjSigma, gLjEpsilon);
//        pForce = &lj;
//    }
//    // start simulation: do job according to job_type
//    std::cout << "Job starts" << std::endl;

//    Dynamics* pDynamics;
//    if (gDynamicsType == DynamicsTypes::kMd)
//    {
//        Md md;
//        md.SetMolSysPtr(mol);
//        md.SetForcePtr(*pForce);
//        pDynamics = &md;
//        // initialize and evolve
//        pDynamics->Initialize();
//        pDynamics->Evolve();
//    }
  
//    // end simulation
//    std::cout << "Job ends" << std::endl;

//    // analysis
//  //  std::cout << "Now computing self diffusion coefficient D ..." << std::endl;
//   // double diff_coef = SelfDiffusionCoefficient("out.xyz", gTimeStep);
//   // std::cout << "D(A^2/fs)= " << diff_coef << std::endl;

//    // @TODO write results to GUI

//    return 0;
//}
